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An  alternative  treatment  of  heat  flow  for  charge  transport  in  semiconductor 
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A  unique  thermodynamic  model  of  Fermi  gases  suitable  for  semiconductor  device  simulation  is 
presented.  Like  other  models,  such  as  drift  diffusion  and  hydrodynamics,  it  employs  moments  of  the 
Boltzmann  transport  equation  derived  using  the  Fermi-Dirac  distribution  function.  However,  unlike 
other  approaches,  it  replaces  the  concept  of  an  electron  thermal  conductivity  with  the  heat  capacity 
of  an  ideal  Fermi  gas  to  determine  heat  flow.  The  model  is  used  to  simulate  a  field-effect  transistor 
and  show  that  the  external  current-voltage  characteristics  are  strong  functions  of  the  state  space 
available  to  the  heated  Fermi  distribution.  ©  2009  American  Institute  of  Physics. 

[doi:  10. 1063/1. 3270404] 


I.  INTRODUCTION 

Throughout  most  of  its  history,  macroscopic  simulation 
of  semiconductor  device  physics  has  focused  on  solving  the 

Boltzmann  transport  equation  using  the  Fermi-Dirac  distri- 

1  2 

bution  function  or  its  Maxwell-Boltzmann  approximation.  ’ 
Representing  mobile  charge  ensembles  with  these  distribu¬ 
tion  functions  is  tantamount  to  treating  them  as  ideal  gases. 
A  three-dimensional  ideal  Fermi  gas  is  spherically  symmetric 
in  momentum  space,  and  its  distribution  in  energy  is  deter¬ 
mined  entirely  by  its  chemical  potential  and  temperature.  In¬ 
serting  the  Fermi-Dirac  distribution  into  the  Boltzmann 
equation  and  solving  for  spatial  variations  in  the  chemical 
potential  and  temperature  form  the  bases  for  drift-diffusion 

.23 

and  hydrodynamic  charge  transport  theories.  ’  These  theo¬ 
ries  offer  accurate  approximations  when  the  rate  that  charges 
interact  within  a  particular  Fermi  gas  is  much  greater  than 
the  rate  that  charges  are  exchanged  between  gases.  An  en¬ 
semble  of  electrons  localized  in  space  and  confined  to  an 
energy  continuum  can  thermalize  very  rapidly  if  the  density 
is  high.  While  devices  often  contain  high  field  regions  where 
the  density  is  low  and  the  distribution  can  become  non- 
Fermi-like,  these  generally  occur  between  high  density  re¬ 
gions  where  the  Fermi  approximation  is  accurate.  If  the  gas 
is  permitted  to  heat  and  the  second  law  of  thermodynamics  is 
enforced  during  the  transition  from  one  Fermi  distribution  to 
another,  the  use  of  Fermi  gases  throughout  provides  a  rea¬ 
sonable  approximation  of  the  electron  dynamics.  Further¬ 
more,  when  the  electron  state  space  includes  forbidden  en¬ 
ergy  gaps  where  no  states  exist,  the  Fermi  gas  approach  is 
flexible  enough  to  accommodate  separate  Fermi  distributions 
for  different  bands  of  states  separated  in  energy.  This  allows 
the  transport  model  to  treat  bipolar  devices  for  which  distri¬ 
butions  of  electrons  in  the  conduction  band  and  holes  in  the 
valence  band  are  represented  by  separate  Fermi  functions 
with  different  chemical  potentials  and  temperatures.  An 
analogous  approach  has  also  been  used  to  accurately  repre¬ 
sent  the  capture  of  injected  electrons  into  the  active  regions 
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of  quantum  well  laser  diodes  by  assigning  separate  Fermi 
distributions  to  the  bound  states  within  the  well  and  the  con¬ 
tinuum  of  propagating  states  above  the  well’s  energy 
barriers.4  The  flexibility,  effectiveness,  and  relative  computa¬ 
tional  efficiency  of  this  Fermi-Dirac  approach  to  charge 
transport  are  evidenced  by  its  widespread  use  in  macroscopic 
semiconductor  device  simulators,  including  computer  aided 
design  software  applications  such  as  PISCES/  DESSIS,6  and 
their  commercial  descendants. 

To  solve  for  spatial  variations  in  chemical  potentials  and 
temperatures,  device  simulators  typically  couple  Poisson’s 
equation  for  the  electrostatic  fields  to  conservation  relations 
for  mobile  charge  densities  and  energies.  Unlike  stochastic 
solutions  of  the  Boltzmann  equation, 7-111  solutions  using  the 
Fermi  distribution  are  only  physically  self-consistent  when 
they  explicitly  enforce  certain  necessary  conditions  of  ideal 
gas  thermodynamics.  One  condition  is  that  heat  flows  be¬ 
tween  gases  with  different  temperatures.  In  a  seminal  paper, 
Stratton1  expressed  this  heat  flow  as  k(T)V  T,  where  k  is  the 
thermal  conductivity  of  electrons  and  T  is  the  electron  tem¬ 
perature.  To  determine  k,  Stratton11  made  use  of  the 
Wiedemann-Franz-Lorentz  law,  an  empirical  observation 
that  relates  the  thermal  and  electrical  conductivities  of  a 
metal  to  its  temperature.  The  great  majority  of  Fermi  gas 
treatments  of  energy  transport  have  adopted  this  concept  of 
an  electron  thermal  conductivity,  and  a  great  deal  of  research 
on  hydrodynamics  has  involved  determining  the  proper 
mathematical  representation  of  this  conductivity  and  the  mo¬ 
ments  of  the  Boltzmann  equation  required  to  complete  the 
model/'1- 

This  paper  proposes  an  alternative  treatment  of  heat  flow 
that  replaces  the  concept  of  an  electron  thermal  conductivity 
with  the  Fermi  distribution’s  heat  capacity.  The  electronic 
equations,  including  electron  density  and  kinetic  energy 
fluxes,  will  be  presented.  Then  it  will  be  shown  that  consid¬ 
ering  the  definition  of  heat  capacity  for  a  Fermi  gas  and 
relating  it  to  these  fluxes  suggests  an  expression  for  heat  flow 
that  completes  the  thermodynamic  model.  To  demonstrate 
this  approach,  a  GaAs  metal-semiconductor  field-effect  tran¬ 
sistor  (MESFET)  will  be  simulated.  Different  parabolic  band 
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structures  will  be  considered  to  show  how  the  interaction  of 
the  Fermi  gases  with  their  state  space  can  affect  the  external 
characteristics  of  the  device. 


13 

ation  time  and  effective  mass  approximations.  '  Neglecting 
any  magnetic  field  and  the  associated  Lorentz  force,  the  first 
moment  for  electrons  in  steady  state  can  be  expressed  as 


II.  DEVICE  EQUATIONS  AND  CHARGE  FLUXES 


A  semiconductor  device  can  be  simulated  by  solving 
equations  for  electrostatics  as  well  as  charge  density  and  en¬ 
ergy  conservation.  Electrostatics  may  be  represented  by  Pois¬ 
son’s  equation 

V  •  eE  =  q{p-n+N+D-N~A),  (1) 


where  e  is  the  dielectric  constant,  E  is  the  electric  field,  q  is 
the  electron  charge,  p  is  the  density  of  mobile  holes,  n  is  the 
electron  density,  N+n  is  the  density  of  ionized  donor  impuri¬ 
ties,  and  N~a  is  the  ionized  acceptor  concentration.  The  carrier 
densities  may  be  computed  by  assuming  Fermi-Dirac  distri¬ 
butions,  constant  effective  mass,  and  parabolic  band  struc¬ 
tures.  For  electrons, 


/o  3/2  roo 

V 2mn 

7 rh3 


:  p 

J  E. 


exp 


(E-F 
knT 


-dE  =  Nc(kBT)i,2Tm(V), 


+  1 


(2) 


where  mn  is  the  effective  mass,  h  is  Planck’s  constant,  E  is 
the  electron  energy,  Ec  is  the  conduction  band  edge,  F  is  the 
chemical  potential  (Fermi  level),  kB  is  the  Boltzmann  con¬ 
stant,  T  is  the  electron  temperature,  and  Nc  is  the  effective 
density  of  states.  The  Fermi  integral  of  order  j  is  defined  as 

ro° 

(3) 


h 


E  ■  Vkf-  W  •  V/ 


dk. 


(7) 


where  the  integral  is  over  all  momentum  vectors  k,  v  is 
electron  velocity,  Tk  is  the  momentum  relaxation  time,  and 
Vkf  denotes  the  gradient  in  momentum  space  of  the  Fermi 
distribution  function  /.  Although  rk  is  generally  a  function  of 
k,  for  simplicity  it  will  be  considered  a  constant.  Also,  a 
simple  constant  effective  mass  band  structure  is  assumed 
such  that  kinetic  energy  E=mn\v\2/2  and  momentum  hk 
=mnv.  The  integral  over  k  can  then  be  transformed  into  an 
integral  over  E,  and  the  flux  can  be  expressed  as 


qTk 

J„  =  -—Nc 

m„ 


E{kBt)mFm{V)  +  (  ) 

3  q 


(8) 


For  numerical  solution,  this  differential  equation  is  repre¬ 
sented  in  discrete  form  on  the  edges  of  a  mesh  that  appro¬ 
priately  represents  the  device.  An  adaptation  of  the 
Scharfetter-Gummel  method1413  can  be  used  to  express  Eq. 
(8)  in  discrete  form.  To  this  end,  Eq.  (8)  is  recast  as  a  differ¬ 
ential  equation  for  electron  density  (2)  as  follows: 


E  +  V 


2Tm  kBT\ 
?>F\n  q  ) 


n 


\  3  /2  q  J 


(9) 


evaluated  for  the  unitless  parameter  77=  (F - Ec)  /  (kBT) .  The 
electron  internal  energy  density  is 


/o  3/2  r°° 

V2  mn 


J  E. 


C E~EC ) 


3/2 


exp(^l  +  1 


dE 


=  Nc(kBT)5l2Jr 3/2(77).  (4) 

Electron  continuity  and  energy  conservation  are  given  by 


If  the  flux  is  assumed  to  be  constant  along  a  given  edge  that 
connects  mesh  points  1  and  2,  then  Eq.  (9)  can  be  solved  for 
the  electron  density  along  the  edge  and  the  corresponding 
electron  flux.  The  resulting  discrete  form  of  the  net  flux  from 
point  1  to  point  2  can  be  written  as 


rl— >2 


qn 

mn 


2^3/2  kBT 
3  T7  !  a  q  . 


-  B(-  £„)n2] 


=  J1(F1,Tl,F2,T2)-J2(Fl,Tl,F2,T2), 


(10) 


dn 

dt 


—  v  •  Jn  +  u, 


(5) 


dEn 

dt 


=  qE-Jn  +  V  ■ST  +  (En)U  + 


En(T)-E„(T,dt ) 
te 


(6) 


where  Jn  is  the  electron  flux  density,  U  is  the  net  electron- 
hole  recombination  rate,  .S'3"  is  the  total  electron  energy  flux, 
( En )  is  the  average  kinetic  energy  per  electron,  7)at  is  the 
lattice  temperature,  and  te  is  the  average  lifetime  for  electron 
energy  loss  from  phonon  emission.  Here  Eqs.  (5)  and  (6) 
have  been  simply  stated  as  conservation  relations,  but  they 
can  be  derived  from  the  zeroth  and  second  moments,  respec- 
tively,  of  the  Boltzmann  equation.  Similar  expressions  exist 
for  mobile  holes. 

Mobile  charge  fluxes  can  be  obtained  from  the  first  mo¬ 
ment  of  the  Boltzmann  equation  using  the  momentum  relax- 


where  L  is  the  length  of  the  edge,  the  notation  (jc)av  denotes 
the  average14  of  quantity  x  evaluated  at  points  1  and  2,  and 
fi(£)  =  £/[exp(£)- 1]  is  the  Bernoulli  function.  For  electrons, 
the  argument  of  the  Bernoulli  function  is  given  by 


2  -T7 3/2  kBT 


2Fy2kBT\ 
3  F \n  q  J 


(11) 


where  £’1_>2  is  the  electric  field  directed  from  point  1  to  point 
2,  and  the  notation  A(jc)  =x2-x1  denotes  the  difference  be¬ 
tween  the  values  of  x  evaluated  at  the  two  points.  Note  that 
net  flux  (10)  appears  as  the  difference  between  electrons 
emitted  from  one  distribution  J  \  and  those  emitted  from  the 
other  J2-  Each  component  is  a  function  of  the  chemical  po¬ 
tentials  and  temperatures  from  both  distributions  because  of 
the  average  prefactor  and  the  arguments  of  the  Bernoulli 
functions.  An  analogous  charge  flux  exists  for  mobile  holes. 
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For  a  Fermi  gas  at  constant  volume,  the  ideal  gas  law 
specifies  the  total  energy  flux  density  ,Sj"1  as  the  sum  of  the 
internal  (kinetic)  energy  flux  density  Skin  and  the  heat  flow 
density  Sheat.  The  internal  energy  flux  can  be  obtained  from 
the  third  moment  of  the  Boltzmann  equation.  To  derive  this 
moment,  the  procedure  outlined  above  is  repeated,  but  the 
Boltzmann  equation  is  multiplied  by  both  electron  velocity  v 
and  kinetic  energy  E  as  follows: 


5  kin 


-ii 


Ev 


Y E  ■  Vkf-  W  ■  V/ 


dk  = 


qni 

m„ 


5£  +  v 

3  \  3  d-3/2  q 


+  I^1ve. 

3  ->=3/2  q 


(12) 


Applying  the  Scharfetter-Gummel  method,  this  flux  is  as¬ 
sumed  constant  along  a  mesh  edge,  reducing  Eq.  (12)  to  a 
differential  equation  for  En  whose  boundary  conditions  are 
the  energy  densities  at  the  edge’s  end  points.  Solving  this 
differential  equation  for  the  flux  gives 
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i-> 

kin 


•2 
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y[B^E)EnA  - 


B(~  &)£„,2)] 
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(13) 


where  Enl  and  Fnl  are  the  energy  densities  at  points  1  and  2, 
respectively,  and  the  argument  of  the  Bernoulli  function  is 
given  by 
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The  heat  capacity  Cel  of  a  Fermi  distribution  deter¬ 
mines  the  energy  £heat  required  to  heat  it  from  one  tempera¬ 
ture  to  another, 


(E  —  F)  ~dk 
dT 


dT 


=  [E,I(T)-Fn(T)]T1fj, 


(15) 


where  T ^  is  the  initial  (final)  temperature.  Similarly,  as 
pointed  out  in  Seeger’s  discussion  of  energy  transport  in 
electron  gases,1  heat  flow  density  is  the  difference  between 
the  internal  and  free  energy  flux  densities.  With  these  ideas 
in  mind,  this  paper  proposes  that  heat  flow  can  be  obtained 
from  Eqs.  (10)  and  (13)  simply  by  evaluating  their  different 
components  at  their  respective  initial  and  final  temperatures 
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(14) 

Analogous  to  electron  flux,  the  net  internal  energy  flux  has 
the  form  of  internal  energy  emitted  from  distribution  1  minus 
that  emitted  from  distribution  2,  and  cross  coupling  occurs 
through  the  prefactor  and  Bernoulli  terms.  A  similar  internal 
energy  flux  can  be  derived  for  mobile  holes. 


III.  HEAT  FLOW  DENSITY 

The  remaining  energy  flux  required  to  complete  the  elec¬ 
tron  thermodynamics  is  the  heat  flow.  Instead  of  invoking  the 
concept  of  an  electron  thermal  conductivity  to  compute  the 
heat  flow,  this  paper  proposes  an  alternative  view.  In  hydro¬ 
dynamics,  charges  moving  through  real  space  are  removed 
from  one  Fermi  gas  with  a  particular  temperature  and  trans¬ 
ferred  to  another  with  a  different  temperature.  When  elec¬ 
trons  move  from  a  Fermi  distribution  at  point  1  to  another  at 
point  2,  they  must  be  heated  or  cooled  from  their  initial  to 
their  final  temperature.  If  T\  <  T2,  electrons  at  point  2  must 
provide  energy  to  the  injected  electrons,  and  this  constitutes 
a  flow  of  heat  from  Fermi  gas  2  to  Fermi  gas  1 .  If  T2<T\, 
the  converse  is  true.  Note  that  the  heat  flow  opposes  the 
temperature  gradient  regardless  of  the  direction  of  electron 
flux. 


By  expressing  heat  flow  as  the  difference  between  inter¬ 
nal  energy  and  free  energy  fluxes  evaluated  at  different  tem¬ 
perature  limits,  including  its  effects  becomes  highly  formu¬ 
laic.  For  any  transport  mechanism  that  can  be  represented  as 
integrals  over  the  energy  spectra  of  two  Fermi  distributions 
exchanging  charges,  the  associated  heat  flow  can  be  readily 
determined  from  Eq.  (16).  For  example,  at  abrupt  heterojunc¬ 
tions  or  Schottky  barriers,  an  ill-defined  electric  field  makes 
drift-diffusion  theory  inappropriate,  and  charge  flux  can  be 
better  represented  by  thermionic  emission.  At  the  conduction 
band  discontinuity,  the  net  electron  flux  from  the  low  L  to  the 
high  H  side  of  the  interface  may  be  represented  as 


J rp  J , 


te  -  J  L-*H  ~  J  H^L  - 


1 

2tt  h 3 


{kBTL)2^ 


1 C,H 


knT, 


-  {kBTH)2Fx 
~Jh(Fh’Th)- 


Fh-Ec,h 

kBTn 


=  Jl(Fl,Tl) 


(17) 


Note  that  the  unitless  arguments  of  both  Fermi  integrals  are 
evaluated  with  respect  to  the  conduction  band  edge  on  the 
high  side  of  the  interface  ECH.  The  corresponding  internal 
energy  flux  is  given  by 
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FIG.  1 .  (Color  online)  Approximation  of  a  GaAs  MESFET.  Dimensions  are 
in  microns.  The  Schottky  barrier  height  is  assumed  to  be  1  eV,  and  each 
contact  is  treated  as  a  perfect  electrical  conductor. 
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Following  the  procedure  indicated  in  Eq.  (16),  the  thermi¬ 
onic  emission  heat  flow  can  be  expressed  as 

Ste,heat  =  SL(FL,  T,  )  -  SL(FL,  TH)  -  (F L  -  Ecm)[Jl(Fl,  Tl) 

-  Jl(Fl ,  Th)]  +  Sh(Fh,  Tl)  -  Sh(Fh,  Th)  -  (Fh 
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(19) 


where  the  chemical  potential  of  each  distribution  is  measured 
from  the  point  of  zero  kinetic  energy,  which  for  thermionic 
emission  is  always  ECH.  The  distribution  on  the  low  side 
feels  the  effects  of  the  barrier  height  through  a  potential  en¬ 
ergy  term  (Ec  H-Ec  L)Jt e,  which  plays  the  same  role  in  its 
energy  balance  equation  as  qE-Jn  plays  in  Eq.  (6). 


FIG.  2.  (Color  online)  Simple  effective  mass  band  structure  to  represent 
GaAs.  Each  band  extends  to  infinite  carrier  energy.  An  effective  mass  of 
0.067ra0  is  used  for  the  T  valley  and  0.85m0  is  used  for  the  L  and  X  valleys. 
The  L  and  X  valley  minima  are  0.29  and  0.48  eV,  respectively,  above  the  T 
valley  minimum. 


lated  by  solving  Eqs.  (1),  (5),  and  (6)  self-consistently  using 
the  full-Newton  method.  Mobile  holes  in  the  valence  band 
and  lattice  heating  were  not  considered.  Fluxes  defined  by 
Eqs.  (17)— (19)  were  applied  at  the  Schottky  gate  interface 
while  fluxes  (10),  (13),  and  (16)  were  used  in  the  GaAs  bulk. 
A  Schottky  barrier  height  of  1  eV  was  used  for  the  gate,  and 
all  contacts  were  treated  as  perfect  electrical  conductors  with 
electron  temperatures  equal  to  300  K.  As  represented  in  Fig. 
2,  conduction  band  electrons  occupied  simple  parabolic  val¬ 
leys  extending  to  infinite  energy.  Different  simulations  in¬ 
cluded  different  valleys.  The  momentum  relaxation  time  rk 
was  treated  as  a  constant,  and  its  value  was  determined  by 
assuming  a  value  of  8500  cm2/  (V  s)  for  the  T  valley’s  low 
held  electron  mobility  //, n=qTklmn  and  assuming  an  effective 
mass  of  m„=0.067m0,  where  m0  is  the  free  electron  mass. 
These  simulations  exhibited  stable  quadratic  convergence  for 
all  applied  bias  conditions. 

For  the  initial  set  of  simulations,  only  the  F  valley  was 
considered.  The  drain  current  versus  drain  voltage  for  differ¬ 
ent  gate  biases  is  shown  in  Fig.  3.  Two  striking  features  of 
these  curves  are  the  superlinear  behavior  at  low  drain  volt¬ 
ages  and  the  low  differential  resistance  at  higher  drain  volt¬ 
ages,  i.e.,  weak  current  saturation  beyond  pinch-off.  Both  of 
these  features  can  be  associated  with  the  constant  momentum 
relaxation  time  Tk  and  the  simplified  band  structure.  At  low 
drain  voltage,  before  the  channel  is  pinched  off,  increasing 
drain  voltage  causes  electrons  to  heat  up,  increasing  their 


IV.  RESULTS  AND  DISCUSSION 

To  test  the  energy  transport  model,  the  simple  GaAs 
MESFET  structure  shown  schematically  in  Fig.  1  was  simu- 


FIG.  3.  (Color  online)  MESFET  drain  current  as  a  function  of  drain  voltage 
for  different  gate  biases.  In  this  simulation,  the  X  and  L  valleys  were  ignored 
and  only  an  infinite  parabolic  F  valley  was  available  to  conduction  band 
electrons. 
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FIG.  4.  (Color  online)  MESFET  drain  current  as  a  function  of  drain  voltage 
for  different  gate  biases.  In  these  simulations,  the  X  and  L  valleys  were 
considered  as  well  as  the  T  valley  in  the  GaAs  conduction  band. 


average  velocity.  Since  rk  is  constant,  Eq.  (7)  shows  that  an 
increasing  average  velocity  combined  with  the  increasing 
electric  field  produces  a  superlinear  drain  current.  At  higher 
drain  voltages,  the  electron  density  in  the  pinch-off  region  is 
reduced,  but  the  high  fields  there  can  heat  the  electrons  to 
well  over  1000  K.  Because  of  the  constant  Tk  and  the  sim¬ 
plified  band  structure,  there  are  unlimited  high  mobility 
states  available  to  the  hot  carriers,  preventing  the  drain  cur¬ 
rent  from  saturating.  In  reality,  rk  should  decrease  with  in¬ 
creasing  electron  energy  and  suppress  both  superlinear  be¬ 
havior  at  low  drain  voltages  and  the  channel  conductance  at 
higher  voltages. 

A  second  set  of  simulations  included  the  X  and  L  valleys 
along  with  the  T  valley  in  Fig.  2.  Although  the  same  constant 
momentum  relaxation  time  rk  was  used,  an  increased  effec¬ 
tive  mass  m„  =  0.85m0  was  assigned  to  the  X  and  L  valleys, 
causing  their  densities  of  states  to  be  higher  and  their  mo¬ 
bilities  to  be  lower  than  the  T  valley.  The  resulting  drain 
current  characteristics  are  shown  in  Fig.  4.  There  is  still  some 
indication  of  superlinearity  at  low  drain  voltages  because 
most  electrons  are  in  the  T  valley  in  this  bias  range,  and  as 
described  above,  the  constant  rk  along  with  increasing  aver¬ 
age  velocity  from  electron  heating  combines  to  produce  this 
effect.  However,  these  simulations  show  increased  differen¬ 
tial  resistance  at  higher  drain  voltages,  and  the  drain  currents 
saturate  at  much  lower  values.  The  saturation  can  be  ex¬ 
plained  by  Fig.  5,  which  shows  electron  density  spectra  at  a 
point  corresponding  to  the  pinch-off  region  of  the  channel. 
The  spectra  were  computed  for  different  drain  voltages  and 


FIG.  5.  (Color  online)  Product  of  the  distribution  function  and  the  density  of 
states  as  a  function  of  electron  energy  at  a  point  in  the  pinch-off  region  of 
the  channel  for  zero  gate  bias  and  different  drain  voltages  computed  with  the 
three  parabolic  valley  conduction  band  model. 


FIG.  6.  (Color  online)  MESFET  drain  current  as  a  function  of  drain  voltage 
for  different  gate  biases  computed  with  both  the  electron  heat  capacity 
model  (dashed)  and  the  electron  thermal  conductivity  model  (solid)  for  heat 
flow.  These  simulations  used  the  infinite  F.  X,  and  L  valley  model  for  the 
GaAs  conduction  band  structure.  The  label  “unstable”  indicates  the  point  at 
which  the  thermal  conductivity  model  diverged  and  produced  no  further 
solutions. 


their  labels  include  the  corresponding  electron  temperatures. 
It  bears  mentioning  that  each  point  in  the  simulation  has  only 
one  electron  chemical  potential  and  temperature,  meaning 
the  T,  X,  and  L  valleys  share  the  same  distribution  function. 
The  different  valleys  only  determine  the  states  available  to 
the  distribution  function.  As  drain  voltage  increases,  heating 
causes  an  increasing  fraction  of  the  total  electron  population 
to  occupy  the  low  mobility  X  and  L  valleys.  Although  the 
distribution  function  is  always  larger  in  the  lower  energy  T 
valley,  the  larger  density  of  states  in  the  X  and  L  valleys 
causes  them  to  contain  the  majority  of  electrons  at  suffi¬ 
ciently  high  temperatures.  Consequently,  as  increasing  elec¬ 
tric  field  in  the  channel  increases  electron  temperature,  more 
carriers  occupy  the  high  energy  X  and  L  valleys,  the  average 
electron  mobility  decreases,  and  the  drain  current  saturates. 

By  way  of  comparison,  the  simulations  that  produced 
Fig.  4  were  repeated  using  the  phenomenological  treatment 
of  heat  flow  Sheat=-/<V T  with  an  electron  thermal  conduc¬ 
tivity  given  by 


k=  q/inn 


T. 


(20) 


The  parameter  v  is  customarily  chosen  to  best  account  for 
the  effect  of  the  electron  scattering  rate’s  energy  dependence 
on  the  Lorenz  number.1 Following  a  similar  simulation 
by  Jyegal  and  DeMassa,”  v  was  set  equal  to  zero.  Note  that 
this  form  for  heat  flow  can  be  applied  only  to  drift  diffusion 
in  the  bulk  regions  of  the  device.  For  thermionic  emission  at 
the  Schottky  interface,  the  temperature  is  discontinuous  and 
its  gradient  ill  defined.  As  a  result,  electrons  on  the  Schottky 
gate  boundary  were  held  at  room  temperature  for  the  electron 
thermal  conductivity  simulations.  Computed  drain  currents 
for  different  drain  and  gate  voltages  are  shown  in  Fig.  6. 
Some  key  features  are  the  more  gradual  onset  of  saturation, 
the  reduced  differential  resistance  after  pinch-off,  and  the 
numerical  instability  exhibited  by  the  electron  thermal  con¬ 
ductivity  model.  For  increasingly  negative  gate  biases,  the 
electron  thermal  conductivity  model  became  unstable  at 
lower  drain  voltages. 

The  drain  currents  in  Fig.  6  can  be  understood  in  terms 
of  the  relationship  between  electron  temperature  and  average 
velocity  that  was  previously  used  to  discuss  Figs.  3  and  4. 
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FIG.  7.  (Color  online)  The  electron  temperature  at  a  point  in  the  MESFET 
pinch-off  region  calculated  for  zero  gate  bias  and  different  drain  voltages 
using  the  infinite  l\  X,  and  L  valley  model  for  band  structure  and  using  both 
the  electron  thermal  conductivity  model  (solid)  and  the  electron  heat  capac¬ 
ity  model  (dashed)  for  heat  flow. 

Figure  7  shows  electron  temperatures  at  a  mesh  point  in  the 
pinch-off  region  computed  for  zero  gate  bias  and  different 
drain  voltages  using  both  the  electron  thermal  conductivity 
and  heat  capacity  models.  The  thermal  conductivity  model 
produces  significantly  more  electron  heating.  At  lower  drain 
voltages,  this  causes  the  lower  mobility  X  and  L  valleys  to 
become  populated  sooner  and  causes  a  reduction  in  the  drain 
current  near  the  onset  of  pinch-off  when  compared  to  the 
heat  capacity  model.  As  indicated  by  Fig.  5,  at  higher  drain 
voltages,  the  electrons  in  the  pinch-off  region  occupy  mostly 
the  X  and  L  valleys.  Therefore,  the  currents  at  higher  drain 
biases  are  dominated  by  the  average  velocities  in  these  val¬ 
leys.  Since  the  thermal  conductivity  model  continues  to  heat 
the  electrons  at  an  accelerated  rate  and  since  the  band  struc¬ 
ture  model  allows  the  X  and  L  valleys  to  extend  to  infinite 
energy,  the  distribution  function  can  access  higher  velocity 
states  and  the  drain  current  continues  to  increase  even  in 
saturation. 

The  differences  between  electron  temperatures  produced 
by  the  electron  thermal  conductivity  and  the  electron  heat 
capacity  models  are  not  altogether  surprising  when  the  ori¬ 
gins  of  the  thermal  conductivity  concept  are  considered.  A 
thermal  conductivity  of  the  form  expressed  in  Eq.  (20)  was 
motivated  by  the  empirical  observation  relating  thermal  and 
electrical  conductivities  in  metals.  A  metal  contains  a  degen¬ 
erate  Fermi  gas  of  conducting  electrons.  From  its  definition 
in  terms  of  its  internal  and  free  energies  (15),  the  heat  capac¬ 
ity  of  a  degenerate  Fermi  gas  is  relatively  low.  However,  the 
MESFET  contains  depletion  regions  under  the  gate  and  in 
the  pinched  off  channel  where  the  electron  densities  are 
many  orders  of  magnitude  less  than  in  a  metal,  and  the  heat 
capacities  per  electron  of  the  rarefied  gases  in  these  regions 
are  much  higher.  Since  the  thermal  conductivity  model  (20) 
contains  no  information  about  the  electron  heat  capacity,  the 
low  electron  density  in  the  depleted  regions  requires  a  larger 
temperature  gradient  to  provide  a  given  amount  of  heat  flow. 
The  alternative  model,  on  the  other  hand,  does  include  the 
enhanced  heat  capacity  of  the  rarefied  electron  gas  and  re¬ 
quires  a  smaller  temperature  gradient  to  achieve  the  same 
heat  flow.  Moreover,  it  might  also  be  argued  that  a  model 
based  on  Fermi  gas  thermodynamics  must  consider  its  heat 
capacity  in  order  to  achieve  self-consistency  and  maintain 
numerical  stability. 


V.  CONCLUSION 

Unlike  stochastic  solutions  of  electron  transport,  solu¬ 
tions  based  on  Fermi  distributions  require  certain  thermody¬ 
namic  conditions  to  be  explicitly  enforced.  One  is  that  heat 
flow  between  Fermi  gases  must  be  included  in  the  total  en¬ 
ergy  flux  in  order  to  fulfill  the  second  law  of  thermodynam¬ 
ics.  For  semiconductor  device  simulation,  computing  this 
heat  flow  has  generally  involved  defining  an  electron  thermal 
conductivity  and  multiplying  it  by  a  gradient  in  electron  tem¬ 
perature.  Determining  the  precise  nature  of  the  thermal  con¬ 
ductivity  has  often  involved  introducing  additional  degrees 
of  freedom  to  the  model  in  the  form  of  higher  order  moments 
of  the  Boltzmann  equation.  This  paper  proposed  a  different 
interpretation  of  heat  flow  between  Fermi  gases. 

Instead  of  invoking  an  electron  thermal  conductivity,  the 
heat  capacity  of  the  Fermi  distribution  was  considered.  This 
can  be  precisely  defined  in  terms  of  the  difference  between 
its  internal  and  its  free  energies.  When  two  Fermi  gases  ex¬ 
change  charges,  the  rates  that  heat  must  be  provided  to  heat 
and  cool  those  charges  can  then  be  determined  by  simply 
comparing  their  kinetic  and  free  energy  flux  densities  evalu¬ 
ated  at  their  initial  and  final  temperatures.  For  energy  trans¬ 
port  based  on  the  drift  and  diffusion  of  Fermi  distributions, 
i.e.,  the  hydrodynamic  model,  this  interpretation  of  heat  flow 
means  moments  0-3  of  the  Boltzmann  equation  are  sufficient 
to  describe  the  charge  thermodynamics. 

In  addition  to  drift-diffusion  transport,  the  same  basic 
scheme  can  be  used  to  determine  the  heat  flow  associated 
with  any  transport  mechanism  that  can  be  expressed  as  inte¬ 
grals  over  the  energy  spectra  of  two  Fermi  gases.  This  was 
shown  for  the  case  of  thermionic  emission  at  an  abrupt  het¬ 
erojunction.  Although  not  implemented  in  this  paper,  the 
method  could  also  be  applied  to  scattering  between  different 
valleys  in  the  semiconductor  band  structure.  If  different  elec¬ 
tron  chemical  potentials  and  temperatures  were  assigned  to 
the  different  valleys,  scattering  between  the  valleys  could  be 
explicitly  considered  by  integrating  the  scattering  rate  over 
the  densities  and  occupation  probabilities  of  the  initial  and 
final  states.  Although  this  is  transport  in  momentum  space 
instead  of  real  space,  it  still  involves  integrating  two  inter¬ 
acting  Fermi  distributions  over  their  energy  spectra.  As  a 
result,  the  same  comparisons  of  internal  and  free  energy 
fluxes  at  the  initial  and  final  temperatures  could  be  used  to 
evaluate  the  heat  exchange  associated  with  the  scattering 
process. 

To  demonstrate  the  model,  a  GaAs  MESFET  structure 
was  simulated  with  simple  parabolic  valleys  representing 
electron  states  in  the  conduction  band.  Each  valley  was  as¬ 
sumed  to  have  the  same  constant  momentum  relaxation  time. 
All  the  simulations  that  used  the  electron  heat  capacity 
model  for  heat  flow  exhibited  stable  quadratic  convergence 
of  the  coupled  device  equations.  They  also  showed  that  the 
terminal  currents  depend  significantly  on  the  states  available 
to  the  heated  Fermi  gases.  This  suggests  that  accurate  band 
structures  and  energy  dependent  momentum  relaxation  rates 
could  be  key  factors  determining  the  accuracy  of  simulations 
based  on  Fermi  distribution  thermodynamics.  Incorporating 
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more  realistic  dispersion  relations  will  be  the  subject  of  con¬ 
tinued  research  into  this  ideal  Fermi  gas  formalism. 
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